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Abstract. We describe a new method for numerical integration of rational 
functions on the real line. Given a rational integrand, we provide a new ra- 
tional function preserving its integral on the line. The coefficients of the new 
function arc explicit polynomials in the original ones. These transformations 
depend on the degree of the input and the desired order of the method. Both 
parameters are arbitrary. The formulas can be precomputed. Iteration yields 
an approximation of the desired integral, with m-th order convergence. Exam- 
ples illustrating the automatic generation of these formulas and a comparison 
with standard numerical schemes are also presented. 



1. Introduction 



The problem of numerical integration of a function over the real line is described 
in Numerical Analysis texts such as [H [TT] . The standard algorithms for the nu- 
merical integration of 



(1.1) 



/ := 



F(x) dx 



start by dealing with the unboundedness of the domain of integration. This is 
usually resolved in two ways: the first one considers the problem on a finite interval 

(1.2) I L := / F{x)dx 



-L 



followed by a convergence study as I — > oo. The alternative is to transform the 
real line to a bounded interval. For example, the map t — x/{\ + x) maps [0,co) 
to [0, 1] and then 



(1.3) 



F{x) dx 



[F(x) + F{-x)) dx 



F 



1 - t 



F 



t - 1 



dt 



{i-ty 



The unboundedness of the original interval of integration is now reflected in the 
(possible) singularity of the new integrand at the boundary t = 1. Observe that if 
the original integrand is a rational function, then so is the new one in (jl.3[) . 

In this paper we present a new numerical method for the integration of rational 
functions on K. It is different in spirit to the standard ones: the numerical approx- 
imation to the integral is obtained from a recurrence acting on the coefficients of 
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the integrand. In particular, the integrand is never evaluated. We illustrate the 
comparison with the standard methods; a more systematic study will be presented 
elsewhere. 

The method presented here is based on a transformation of the coefficients of a 
rational function, that preserves its integral. This is the so-called rational Landen 
transformation. The original problem discussed by Landen, Gauss, Legendre and 
others deals with the elliptic integral 

^l 2 al9 

(1.4) G(a,b)= 1 



V a 2 cos 2 9 + b 2 sin 2 6 

Gauss [5] made the numerical observation that the function G(a, b) was invariant 
under the transformation 

(1.5) (a,*)- 



2 ' 

The iteration of (|1.5|) leads to a numerical evaluation of the elliptic integral, or, as it 
has been explained by J. and P. Borwein in [3] , to efficient methods for the numerical 
evaluation of tt. The rational analogue of this transformation was developed in 
[H UOj and here we show how to use it as a numerical method to evaluate rational 
integrals. The reader will find in [9] a survey of the diverse aspects related to these 
transformations. 

Section 2 discusses the basic structure of the algorithm. Section 3 introduces 
a family of polynomials that play an important role in the development of the 
formulas given in Section 4. Finally, Section 5 contains some examples. The first 
illustrates the steps for a method of order 2 acting on a rational function of degree 
6. The next two examples illustrate the accuracy of the method and its comparison 
to the trapezoidal rule. A systematic study of the cost involved in this algorithm 
will be presented elsewhere. 

2. The Landen transformation and algorithm 
We present a general description of an iterative algorithm for the evaluation of 

POO 

(2.1) /:= / F(x)dx. 



Here F is a rational function given as 

(2-2) F(x) = §4' 

A{x) 

with 

p p-2 



(2.3) A(x) = a kX p ~ k and B(x) = ^ b k 

k=0 k=0 

where a p =/= 0. The general construction treats the coefficients aj and bj as in- 
determinates. Naturally, for specific integrands the parameters a, and bj are real 
numbers and the maps described in this section are defined on parts of M. 2p where 
the integrals are convergent. 

The set 

lZ 2p : {(do, ai, . . . , a p , bo, ... , b p -2) '■ such that F in (|2.2[) has finite integral }, 



NUMERICAL LANDEN 



3 



will be used to represent the rational function F in terms of its coefficients. It will 
be referred as the coefficient space. 

In Section d] we describe the construction, for each integer m > 2, of a rational 
function Fi m that satisfies 

/oo />oo 
F hm (x)dx = / F{x)dx. 
-oo J — oo 

The function F\^ m {x) has the same degree as the original F and the new coeffi- 
cients are polynomials in the old coefficients <2q, . . . , a p , bo, ■ ■ ■ , 6 p -2- Naturally, this 
produces a map on lZ 2p that we denote by £ m ,p '■ 7Z 2p — > Ti 2p , called the rational 
Landen transformation of order m and degree p. We denote by £^ „ its n-fold com- 
position. Now introduce a new operator by : TZ 2p — > K as the map which takes a 
vector of length 2p corresponding to a rational function in 1Z of degree p and returns 
the value of the rational function at x = 0. In terms of the vector of coefficients, 
this is the ratio of the last entry over its (p + I)— th one. The composition <j> o £^ p 
will be denoted by 0^ p . The motivation behind <fi is the following: if 

(2.5) F[x) _b + b 1 x + --- + b p - 2 x p - 2 



oq + a\x + • • • + a p x p 
then the function obtained by iterating the map £ miP n times is written as 

(2.6) F n , m {x) = ^ + K» X + --- + bp -^ xP ;\ 

ao, n + ax, n x H h a p ^ n x p 

The value of F n ^ m at x — gives a sequence of real numbers that converges to 1/tt 
times the integral of F in (|2.5[) : see [10] for details. The update on this sequence 
comes from the coefficients of F n ^ m . These are obtained by applying £ m . p to those 
of F n - lt7n . Finally, define a := (a , a p , b , 6 p _ 2 ). 

The main result of [10] is that the Landen transformation satisfies 

(2-7) C, P («) ^ - 

as n — > oo if / < oo. Furthermore, the convergence is of order m, that is, 

j m 



(2.8) 



< C 



C,p(«) 



This convergence result appears in [B] for the case m = 2 and the general case can 
be established along the same lines. See [5] for details. 

Adapting the rational Landen transformations into a numerical method for cal- 
culating / involves a process with two parts. The first one is a symbolic calculation 
of the explicit algebraic formulae for the rational Landen transformation. The 
steps in this calculation are described in detail in Section @] The second one is the 
iteration of these formulas. 



Algorithm 1 
Input: 

1) An integer p: the degree of the denominator A. 

2) An integer m > 2: the order of the transformation. 

Output: 
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The explicit formula for the transformation £ m , P ■ r R? p —* lZ 2p , as polynomials 
in the a, and bj. 

Note 2.1. Observe that, given m and p, the map £ mjP can be precomputed and 
the result can be stored for its use in the second algorithm. Therefore, the first 
algorithm carries a one-time cost and is not figured into the time of the method. 
This precomputation will be assumed in the discussion of the second algorithm. 

Algorithm 2 
Input: 

1) A vector a representing the coefficients of the rational integrand F. 

2) An integer m > 2: the order of convergence. 

3) An integer n G N: the number of iterations of the Landen map. 

Output: The expression <fi™ lp (d) that approximates I/tt. 

3. The evaluation of the polynomials P m and Q m 
The algorithm described in Section [4] employs the polynomials 

|m/2J . . 

(3-1) P m (x):= J2 (-^{^y™-* 

and 

L(m-1)/2J 

(3-2) Q m {x):= J2 (- 1 )'( 2j + 1 )^ (2i+1) - 

The integration algorithm discussed here is based on the fact that the rational 
function 

(3.3) R m ( X ) := ^ 
satisfies 

(3.4) cot {me) = Rm (cot9). 
For instance, for m = 2, we have 

(3.5) P 2 {x) = x 2 -l, and Q 2 (x) = 2x. 

Note 3.1. This trigonometric property is instrumental in the proof that the integral 
of F is the same as that of F\. See [TU] for details. 

4. The algorithm 

In this section we describe each of the steps in the first algorithm. This algorithm 
has been implemented in Mathematica 6.0. 

Step 1. Construct the polynomial 

p 

(4.1) H(x) =J2 h * xP ~ 1 

i=0 

defined by 

(4.2) H{x) := Res, (A(z), P m (z) - xQ m (z)) , 
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where Res z denotes the resultant in the variable z. The degrees of the polynomials 
involved in ()4.2j) are p = degA and m = deg(P m (z) — xQ m (z)), respectively. 

The resultant. Given two polynomials a(t) and f3(t), the resultant of a and [3 is 
defined by 

r s 

(4.3) Res(a,/3) := JJ JJCy^ - a;,-), 

i=ij=i 

where a;, are the roots of a(t) = and are the roots of j3{t) = 0. 

The resultant of two polynomials can be computed as the determinant of the 
Sylvester matrix formed by their coefficients; see [7]. For instance, if 

a(t) = a + ait + a 2 t 2 + a 3 t 3 and (3(t) = b + bit + b 2 t 2 , 

then the Sylvester matrix is defined by 



/a3 


a 2 


«1 


a 










as 


(12 


«i 


«o 













(12 


«i 


a Q 


b 2 


h 


b 














&2 


h 


b 














&2 


h 


&o 













&2 


h 





and it is a square matrix of size deg(a) + deg(/3) + 2 = 7. The resultant of a(t) and 
(3{t) in ga|| is 

Res(a, /3) = — a 2 a3&o6i + aia 3 6 6i — ao a 3^i + a 2^o^2 

-2aia a blb 2 - aia 2 6o^i92 + 3a a 3 6o6i62 + a a2blb 2 
+ai&o&2 — 2aoa2bo&2 — 00016162 + 0^2 ■ 
In general, the resultant of two polynomials is a polynomial in their coefficients. 

Note 4.1. The polynomial H has the same degree as A, the denominator of the 
integrand R(x). It will become the denominator of the new rational function. Its 
coefficients hi are polynomials in those of A. The calculation of H can be obtained 
by evaluating the determinant of a square matrix of dimension deg(A) + deg(P m ) + 
2 = p + r7i + 2. For instance, for m = 2 and 

(4.5) A{x) = aox* + aix 3 + a2X 2 + a\X + 04, 
we obtain 

H{x) = 16aoa4X 4 + 8(ai&4 — aoa3)a; 3 + 4(aoa2 — 0103 + 4aoa4 + a204)a; 2 

+ 2(— &oai + &i<22 — 3aoa3 — 02(13 + 3ai<24 + a 3 a^)x 

+ (a — di + a 2 — a 3 + a 4 )(a + ai + a 2 + a 3 + a 4 ). 

Step 2. Form the polynomial 

(4.6) E( X ) =Y, h i {Pm{x)Y~ l (QMY . 

i=0 
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Step 3. The polynomial E{x) formed in Step 2 is a multiple of the denominator 
A(x). Compute the quotient 

m »(,) - 1| 

and write it as 

r 

(4.8) Z(x) = Zkx r ~ k , with r=p(m- 1). 

Step 4. Compute the product 

(4.9) C(x) = B(x)Z(x) 
and write it as 

(4.10) C(x) = c k x s - k , with s = mp- 2. 
Step 5. Form the expression 

(4.11) TM,I,U fj L)° ' '(.,.",)(')• 

for a, b, x e N. 

Step 6. Define the expressions 



x [T x+am (2j,s-2j) + (2i,s-2j)], 



2 2(a-/3) a - /V - a - 1 + 

7 



and 



M 2 (j,a,/3, 7 ,m,f>) := (-l)^c 2j+1 2 2 ^ * 

X [T A+am (2j + 1,5- 2j - 1) - Tx—am (2.7 + 1, * - 2j - 1)] , 



with i/ := p/2 and A := (mp — 2)/2. 
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Step 7. Define 



v—X / ,\ a 



7=0 \ V 1 7 J=0 

!/— 2 / A 1 — 7 a \ 

7 __0 \ j=0 a=l /3=0 / 
^ v — 1 / A v— 1 a \ 

+ yE E E E M 1 (j ! a,A7,m,p) 

7=1 \ j=0 a=y— 7 /3=a — J/+7+1 / 
!/-2 / \-X v-X-y a-X \ 

+ ?E EEEw-^r 1 

7 =0 W=0 a=l /3=0 / 
i/-2 /a-1 w— 1 a—X \ 

+ iE EE E M 2(j,a,/?, 7 ,m,p) 



2 s 

7=1 \ j=0 a=J/— 7 0=0 



Step 8. The new rational function is denned by 

J(x) 



(4.12) F hm {x) 



H(x) 
It satisfies (j2~l)l . 

The reader is referred to [TU] for the proofs of the formulas describing this algo- 
rithm. 



5. Examples 

In this section we give examples that illustrate the rational Landen transforma- 
tions. 

Example 5.1. We provide a step by step construction of the rational Landen 
transformation of order m = 2 for the function 

b x 4 + hx 3 + b 2 x 2 + b 3 x + bi 

(5.1) F(x) 



aox 6 + ciix 5 + a 2 x 4 + a^x 3 + a^x 2 + a§x + 
The goal is to produce a function Fi^(x) with the same integral as F(x). 



Note 5.1. The special case &i = b% = a\ = a 3 = — was the first example of 
this new type of transformation. It appears in pQ. 

Note 5.2. The choice of m = 2 requires the evaluation of the polynomials P2(x) — 
x 2 — 1 and Q 2 (x) = 2x. These are computed directly from (|3.1[) and 
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Step 1 computes the polynomial H(x) from (|4.2p . The Mathematica command 
Resultant yields the expressions 



ho = 64aoae 

hi = — 32(ao<25 — aide) 

hi — 16(aoa4 — 0105 + 6aoa6 + a2i6) 

hs = — 8(aoa3 — 0104 + 5aoa5 + 0205 — 5aia 6 — a 3 a 6 ) 

hi = 4(aoa2 — 0103 + 4aoa4 + 0204 — 4ai<25 — 0305 + 9oqOq + 4a2a6 + 0,4(16) 

/15 = — 2(aoai — aia2 + 3ooa3 + 0203 — 3etiet4 — 0,304 + 50005) 

— 2(3a2<25 + 0405 — 5etiQ,6 — 30305 — a^ag) 

h e = (a — a>\ + a 2 — 03 + 04 — a 5 + a e )(a + a\ + a-i + a 3 + a 4 + 05 + a 6 ). 



The polynomial H(x) is the denominator of the new rational function obtained as 
a product of the Landen transformation. 

Step 2 computes the polynomial E{x) from (|4.6p . In this example, this is a poly- 
nomial of degree 12 (= mp), that we write as 

12 

(5.2) E{x) =J2 e ^ 12 ~ l - 

i=0 

The symbolic expansion of (|4.6[) produces 

e = e 12 = 64a a 6 

ei = -en = 64(a a 5 - a x a 6 ) 

C2 = eio = 64(aoa 4 — 0405 + a 2 ae) 

e3 = -eg = 64(a a 3 - aia 4 + a 2 a 5 - 030 6 ) 

64 = es — 64(aoa 2 — aia.3 + 0204 — 0305 + 0.40.0) 

e§ = — e 7 = 64(aoOi — 0402 + a 2 03 — a 3fl4 + 0405 — 05^6) 

r* a ( _ 2 2i 2 2 1 2 2i 2 \ 

e6 = o4(a — a x + a 2 — a 3 + a 4 — a 5 + a 6 ). 

Step 3 computes the quotient of E(x), produced in Step 2, and A{x), the denom- 
inator of the original integrand. In the example discussed now we obtain 

(5.3) Z{x) = 64(a e x 6 — a$x 5 + a^x 4 — a 3 x 3 + a 2 x 2 — a\x + a ). 

Step 4 simply evaluates the product C(x) — B(x)Z(x), where B(x) is the numer- 
ator of the original integrand and Z{x) comes from Step 3. The polynomial C{x) 
is written as 

10 

(5.4) C(x)=Y, c * xl °~ k > 

fc=0 
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with 



CO 


= 64a 6 6 






Cl 


= -64(a 5 6 


- a 6 6i) 




C2 


= 64(0460 — 


a 5 h + a 6 b 2 ) 




C3 


= -64(a 3 6 


- a 4 bi + a 5 b 2 


- a 6 b 3 ) 


C 4 


= 64(a 2 6 — 


0361 + a 4 b 2 - 


a 5 b 3 + a 6 b 4 ) 


C5 


= -64(ai6 


— a 2 b\ + 0362 


- a 4 b 3 + 0564) 


C6 


= 64(a 6 - 


a\b\ + a 2 &2 — 


a 3 b 3 + a 4 b 4 ) 


C7 


= 64(a foi - 


aife 2 + a 2 b 3 - 


a 3 b 4 ) 


C8 


= 64(a 6 2 - 


a\b 3 + a 2 b 4 ) 




eg 


= 64(a 6 3 - 






ClO 


= 64a 64. 







Step 7 combines the functions defined in Step 5 and 6 to produce the new numer- 
ator 

4 

(5.5) JW = ^4^", 

fc=0 

with 

d a = 32(a 6 6 + a b 4 ) 

d 4 = —16(a 5 b — cisbx + a b 3 — a 4 b 4 ) 

d 2 = 8(a 4 b + 3a 6 b — a^bx + a b 2 + a e b 2 — + 3a b 4 + a 2 b 4 ) 
d 3 = -4(a 3 6 + 2a 5 6 + a &i — a 4 bi — 2a 6 b 4 — a\b 2 + a 5 6 2 ) — 

-4(2a 6 3 + a 2 b 3 - a 6 b 3 - 2a 4 b 4 - a 3 b 4 ) 
d 4 = 2(a 6 + a 2 b + a 4 b + a 6 b - a 4 bi - a 3 b 4 - a 5 b 4 ) + 

+2(a & 2 + a 2 b 2 + a 4 b 2 + a 6 b 2 - a 4 b 3 - a 3 b 3 + a b 4 + a 2 b 4 + a 4 b 4 + a 6 b 4 ). 



Note 5.3. Given a rational function of order p and a choice of method of order 
to, the calculation of H and J illustrated here is done once. We have produced a 
transformation sending 

(5.6) F(x) = b ° x4 + + ^ + hX + h 

a x e + aix 5 + a 2 x 4 + a 3 x 3 + a 4 x 2 + a 5 x + a e 

to 

d^x A + d\x 3 + d 2 x 2 + d 3 x + d 4 



(5.7) Fi >6 (x) := £ 6 , 2 F(x) 



hgx 6 + h 4 x 5 + h 2 x 4 + h 3 x 3 + h 4 x 2 + /i 5 x • 



with the new coefficients given as above. We expect to produce a precomputed array 
of formulas, indexed by (p, to), to be made available to the community. 



In the next series of examples, we will assume that the formulas for the Landen 
transformations have been precomputed. 
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Example 5.2. The rational function 

(5.8) F (aO = 
satisfies 

r°° tv 

(5.9) I := / F (x)dx 



11 

with numerical value 

(5.10) I - 0.94722582509948293643. 

We now employ the algorithm, with a method of order 2, to obtain the sequence of 
rational functions 



(5.12) /= / F n>2 (x)dx 



(5.11) F n , 2 (x) = - , . 

x £ + a n x + b n 

with Fo t 2(x) — Fq(x), and the property that 

dx 

The convergence analysis described in [TU] shows that a n — > 0, b n — * 1, thus 
(5.13) I = lim 7r c„. 

n — ^oo 

Even though the limiting value of the integral depends only upon the terms c„, the 
formulas to generate these values also involve a n and b n . Therefore one must store 
the current value of all the parameters. The first few of them are shown in the next 
table: 



n 




a n 


bn 





l 


4 


15 


1 


8 


28 


4 


15 


15 




2 


1 


7 


4841 




3 


10 


3600 


3 


8441 


8687 


64900081 




29046 


96820 


69710400 



Table 1 . Rational Landen of order 2 



Note 5.4. The expression 7rc„ gives an approximation to the integral of F (x). For 
example, for n = 6, we find that 

3471070386673821384824326347489289738211683509253931254760471 
° 6 ~ 11512238093504492278949475398059063785494372327433955614454608 
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and ttcq agrees with the integral of Fq up to 13 digits. The coefficients (a n , b n , c n ) 
are rational numbers and their height increases with n. Recall that the height of a 
rational number, written in irreducible form as x — — , is 

(5.14) h(x) = Max{|u|, \v\}. 

For example, at n = 10, the heights of a%o, bio, cio have approximately 1000 digits. 
At this stage, the value ttcio gives 196 digits of the integral. Naturally, this increases 
the complexity of the calculations if we use exact arithmetic. An interesting way 
out of this problem is to replace the rational number c„ by the truncation of its 
continued fraction. For example, the first 20 terms of the continued fraction of cio, 
a rational number of height only 6, differs from c%o by less than 10~ 24 . Details will 
be given in a future publication. 

Example 5.3. Mathematica 6.0 shows that if 

(5.15) F(x) ' 
then 



(5.16) I := [ F(x)dx = 2TrJ-?-{V37-5). 

J —oo V 1 1 1 

We employ the identity (|1.3[) to map the problem to the interval [0, 1]. The new 
rational function 

60a; 6 - 288a; 5 + 584a; 4 - 648a; 3 + 422a; 2 - 156a; + 26 
9 ^ X ' '~ (3a; 4 - 15a; 3 + 31a; 2 - 31a; + 13) (57a; 4 - 153a; 3 + 157x 2 - 73a; + 13) 

60a; 6 - 288a; 5 + 584a; 4 - 648a; 3 + 422a; 2 - 156a; + 26 
171a; 8 - 1314a; 7 + 4533a; 6 - 9084a; 5 + 11485a; 4 - 9314a; 3 + 4707a; 2 - 1352a; + 169 

satisfies 

/>1 />oo 

(5.17) / g(x)dx = / F(x)dx. 

JO J -oo 

We now compare the numerical approximation to the integral of g over [0, 1] com- 
puted with the methods described here, with a numerical integration using the 
trapezoidal rule. A more systematic comparison with more sophisticated classical 
numerical schemes will be described elsewhere. 
The trapezoidal rule states that 

(5-18) J" g(x) dx = ^(g(a) + g(b)) - ^V'(0, 

where h = b — a and £ G [a, b]. Define 

(5.19) M := Max{|s"(t)| : a < t < b}. 

Then the error term, namely the second term in (|5.18[) is bounded by Mh 3 / '12. To 
obtain an approximation to the integral I in (|5.16[) . we choose n G N, partition 
[0, 1] into n intervals of equal length h = 1/n and apply the trapezoidal rule to 
each subinterval. This yields the expression 



(5.20) /„ := 



„1 / n-2 \ ^ n— 1 

l 9{x) dx ^^[ 9o + T, 9i + 9n ) - 12^ E 
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where g L — g(i/n) and i — = < £j < — . The total error in this approximation is 
bounded by M/12n 2 . To compute the approximation I n to / requieres the n values 
{gi ■ < i < n}. The relative error (I - I n )/I for n = 100 is 5.29805 x 10~ 6 . It 
drops to 3.1505 x 10~ 8 for n = 1000 and to 2.9445 x 10" 10 for n = 10000. 

We now use the method of rational Landen transformations of order m to produce 
approximations to the integral / of F(x) over R. Recall that the method yields a 
family of rational functions R n , m {x) with integral /. For example, the first two 
functions for a method of order 2 are 

4(2a; 2 + 6a; + 15) 



R 1>2 (x) = 

i?2, 2 (x) = 



208x 4 + 456a; 3 + 600a; 2 + 396a; + 171 ' 

8(13848a; 2 + 11652a; + 11531) 



569088a; 4 - 35136a; 3 + 756384a; 2 - 8616a; + 232537' 
The approximations to I are then obtained from 

„„ . Constant term in the numerator of R n m (x) 

n,m Constant term in the denominator of R n .m(x) 

The next table shows the relative errors 

(5-22) r eln. m := jyj 

for 2 < m < 6 and 2 < n < 5. 



n 


m = 2 


m — 3 


m = 4 


m = 5 


m = 6 


2 


0.30314 


0.022076 


0.0021170 


2.2646 x 10~ 6 


6.3257 x 10~ 7 


3 


0.058475 


0.000035272 


5.2932 x 10" 12 


2.9440 x 10" 23 


4.4813 x 10~ 40 


4 


0.0021170 


3.2713 x 10" 15 


2.0616 x 10" 47 


1.9758 x 10" 115 


3.6655 x 10" 239 


5 


3.2700 x 10" 6 


3.6952 x 10" 45 


5.3750 x 10" 190 


3.1671 x 10" 577 


4.0442 x 10" 1434 



Table 2. Relative error for the numerical evaluation of /, 



This table contains clear evidence to support the convergence orders claimed in 
(EH). 



Example 5.4. The rational Landen transformations can be used to evaluate 

f°° dx - n 
(5 - 23) J-oo {x-2y + e>-l> 

for e > small. This example illustrates the fact that the proposed method con- 
verges, even when the integrand has poles very close to the real axis. A systematic 
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description of the sensitivity of the iteration with respect to the parameter e, will 
be presented elsewhere. 

For fixed e > 0, we apply a method of order 2 to (|5.23p . This produces a sequence 
of rational functions of the form 

(5.24) R n (x) 



that satisfy 

/CO 
R n (x)dx = -. 
-oo e 

The explicit Landen transformation of order 2 is given by 

(5.26) ao,n+l = (ao,n - a l,n + Ct2,n){ao,n + cti jn + a 2 , n ) 

Ol,n+l = 2ai !n (ao,n — Cl2,n) 

&2,n+l = 4a 

^0,n+l = 2&o in (ao, n + a 2,n)i 

with initial conditions 

(5.27) a , = 4 + e 2 , ai )0 = -4, a 2 , = 1, &o,o = 1. 

The sequence R n (x) has coefficients that depend upon the parameter e. For 
example, 

D, _ 2(5+e 2 ) 

"IW — 4(4+e 2 ):z 2 -8(3+£ 2 ):r + (l+e 2 )(9+e 2 ) 

and 

p / n, _ 4(5+e 2 )(25+14e 2 +e 4 ) 

-n-2 W — (l+6e 2 +£ 4 )(49+22e 2 +£ 4 )-16(-l+e)(l+e)(3+e 2 )(7+e 2 )x+16(l+e 2 )(4+c 2 )(9+e 2 ) ■ 

The theory described above shows that, for fixed e > and n — > oo, the sequences 

(5.28) ^^L, ^-0, 

U0,n Oo.n Oo.n 

converge to the stated limits. Moreover, the invariance of (|5.23[) under the trans- 
formations given in (|5.26p show that L = e. 
Define the error 

2\ V2 



(5-29) err„ := _ e ) + (:_) + 



Then Table 3 shows the ratios erri6/erri5 obtained after 15 iterations of (15.26|) for 
methods of order 2 and 3. The calculations are done with 10 6 digit precision. 

The data in Table 3 shows the exponential decay of the error. Given a tolerance 
5 > 0, we have observed that the number of steps required to achieve err„ < 5 
increases as e — > 0. A quantitative description of this phenomena is in preparation 
and it will be reported elsewhere. 
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UI ClCl z 




.1 


3.58047 x 10" 569 


3.49118 x 10 - 497802 


.01 


1.36862 x 10" 57 


4.24935 x io~ 49853 


.001 


2.07254 x 10~ 6 


4.73905 x 10" 4986 


.0001 


2.16805 x 10~ 2 


3.48094 x 1CT 499 


.00001 


4.68150 x 10" 1 


1.62880 x 10~ 50 



Table 3. The quotient errig/erris as a function of the parameter e. 



6. Conclusions 

We have described the rational Landen transformations and their use in the 
numerical integration of rational functions. We have exhibited fast convergence of 
this method and presented an example comparing it to the classical integration 
schemes. 

A systematic comparative analysis of this method with respect to standard nu- 
merical algorithms will be discussed elsewhere. An interesting challenging problem 
is to extend the use of rational Landen transformations to produce fast numerical 
integrators for arbitrary functions. In particular the method is well suited for the 
numerical integration of meromorphic function with poles off the real line. 
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